Mon. Not. R. Astron. Soc. 000, 1-21 (2012) Printed 3 March 2013 (MN WT^ style file v2.2) 



A dynamical analysis of the Kepler- 11 planetary system 



Cezary Migaszewski , Mariusz Slonina t and Krzysztof Gozdziewski^'^$ 

^Torun Centre for Astronomy, Nicolaus Copernicus University, Gagarin Str. 11, 87-100 Torun, Poland 



Accepted 2012 August 23. Received 2012 August 20; in original form 2012 May 3 



ABSTRACT 

The Kepler- 11 star hosts at least six transiting super-Earth planets detected through the pre- 
cise photometric observations of the Kepler mission (Lissauer et al.). In this paper, we re- 
analyze the available Kepler data, using the direct A/^-body approach rather than an indirect 
TTV method in the discovery paper. The orbital modeling in the realm of the direct approach 
relies on the whole data set, not only on the mid-transits times. Most of the results in the orig- 
inal paper are confirmed and extended. We constrained the mass of the outermost planet g to 
less than 30 Earth masses. The mutual inclinations between orbits b and c as well as between 
orbits d and e are determined with a good precision, in the range of [1,5] degrees. Having sev- 
eral solutions to four qualitative orbital models of the Kepler- 1 1 system, we analyze its global 
dynamics with the help of dynamical maps. They reveal a sophisticated structure of the phase 
space, with narrow regions of regular motion. The dynamics are governed by a dense net of 
three- and four-body mean motion resonances, forming the Arnold web. Overlapping of these 
resonances is a main source of instability. We found that the Kepler- 11 system may be long- 
term stable only in particular multiple resonant configurations with small relative inclinations. 
The mass-radius data derived for all companions reveal a clear anti-correlation between the 
mean density of the planets with their distance from the star. This may reflect the formation 
and early evolution history of the system. 
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1 INTRODUCTION 

The Kepler space mission is a breakthrough in the field of searches 
for the Earth-like extrasolar planets (Borucki et al. 2010; Koch 
et al. 2010; Jenkins et al. 2010; Caldwell et al. 2010). About of 
150,000 solar dwarfs are monitored by 0.95 — meter Kepler tele- 
scope. The photometric data are publicly available from the MAST 
archive^ . 

To date, the mission identified more that 2,200 planetary can- 
didates (Batalha et al. 2012). Among them, many multi-planet sys- 
tems are found. For instance, planets were confirmed in two-planet 
configurations, i.e., Kepler-10 (Batalha et al. 2011; Fressin et al. 
201 1), Kepler-25, 26, 27, 28 (Steffen et al. 2012), Kepler-29, 31, 32 
(Fabrycky et al. 2012), Kepler-23, 24 (Ford et al. 2012); in three- 
planet systems Kepler-9 (Holman et al. 2010), Kepler-30 (Fabrycky 
et al. 2012), Kepler-18 (Cochran et al. 2011); in four-planet sys- 
tems (Borucki et al. 2011), as well as in five-planet configurations 
Kepler-20 (Gautier et al. 201 1 ; Fressin et al. 201 1), Kepler-33 (Lis- 
sauer et al. 2012). The Kepler- 11 hosts six planetary companions 
(Lissauer et al. 201 1). The transiting planet candidates can be con- 
firmed through determining their masses with the help of the so- 
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called Transit Timing Variations method (TTV, Holman & Mur- 
ray 2005; Agol et al. 2005). In this approach, the (O-C) variations 
between observed mid-transit times and their ephemeris are the 
observables, which can be fitted by an appropriate orbital model. 
In recent papers, also additional observables are analysed, like the 
so-called Transits Duration Variations (TDVs) (see, e.g., Nesvorny 
et al. 2012). 

In this paper, we re-analyse the photometric data of Kepler- 1 1 
with a modified, direct approach providing an alternate estimation 
of masses and orbital elements. To describe this method further 
in the paper, we recall shortly the main conclusions in (Lissauer 
et al. 2011). Using the TTV method and an assumption of strictly 
coplanar model of the system, they determined masses of five inner 
planets in the range of a few Earth masses. The outermost planet 
interacts weakly with the inner companions, and its mass could be 
roughly constrained as smaller than the Jupiter mass. It has been 
not confirmed as a planet, although the probability of blending is 
very small, ^0.001. Orbital eccentricities in the Kepler- 11 system 
were determined only for the five inner objects. Due to the assump- 
tion of coplanarity, a determination of mutual inclinations between 
the orbits was not possible. Lissauer et al. (2011) argue that these 
inclinations should remain in the range of [0,2] degrees. The dy- 
namical analysis have revealed that the system is not involved in 
the mean motion resonances (MMRs), however a pair of planet b 
and planet c is close to 5:4 MMR. 

The determined masses and radii of the planets imply con- 
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strains on their chemical composition. Planets d, e and f might have 
similar internal compositions to those of Uranus or Neptune, while 
planets b and c are rather ice rich, with a smaller amount of H2/He 
mixture than these planets in the Solar system. 

In this paper, we focus mostly on the global dynamics of the 
system and a few aspects which were not addressed in the discovery 
paper. 

First of all, we model the available Kepler data through a di- 
rect algorithm that relies on the self-consistent A/^-body fitting of the 
light-curves, instead of the TTV method applied in the discovery 
work. The TTV algorithm makes use of the transit times a poste- 
riori, after they are determined from the light curves. Through ex- 
tensive numerical experiments, we found that the direct approach 
brings more information than the TTV method. For instance, we 
could constrain the mass of the outermost planet to less than ~ 
30 Earth masses. We also found significant bounds for the mutual 
inclinations to less than 5° for planets b and c as well as for plan- 
ets d and e. 

The direct model, also called the dynamical-photometric 
model, already was used in a few papers. For instance, it was ap- 
plied to analyse the light curve of the triple-star system KOI- 126 
(Carter et al. 2011), and to estimate masses of two planets transit- 
ing Kepler-36 (Carter et al. 2012). This algorithm also verified the 
Kepler-9 model, which was found first with the help of the TTV 
algorithm (Holman et al. 2010). 

A number of initial conditions found with the direct approach 
makes it possible to investigate the dynamics of the system. We 
focus on the short-term time scale, governed by the mean mo- 
tion resonances. We study the multi-dimensional structure of the 
phase space with the help of dynamical maps. In the vicinity a few 
qualitative transit models considered in this work, the dynamics are 
governed by a dense net of 3-body and 4-body mean motion reso- 
nances. This net may be identified with the Arnold web, which is a 
feature of close to integrable Hamiltonian systems. The Kepler- 11 
appears as strongly resonant extrasolar system, and this feature may 
reflect its trapping into MMRs at the early stages of the formation 
and evolution. 

Using a new determination of the masses and radii, we found 
a curious mass-radius relation implying a clear anti-correlation be- 
tween the mean density of the planets and their distances from 
the star. Their densities exhibit a sequence of planet b which is 
denser than Neptune, through the Neptune-like planet c. Uranus- 
like planet d, Jupiter-like planets e and f, and planet g which is 
likely Jupiter/Satum-like. 

The paper is structured as follows. In Sect. 2, we shortly de- 
scribe the photometric data of Kepler- 11 available in the MAST 
archive. We also refine the observational TTV model. In Sect. 3, we 
present the results derived through intensive computations with the 
bootstrap algorithm. Furthermore, we discuss a possible composi- 
tion of the planets (Sect. 4). Section 5 is devoted to the dynamical 
analysis of the Kepler- 11 system. Conclusions and prospects for a 
future work are given at the end of this paper. 



2 TRANSITS IN A MULTI-PLANET SYSTEM 

The photometric data of Kepler- 11 were taken from the MAST 
archive. At the time of writing this paper, the publicly available 
light-curves span about of 500 days in six parts. These data were 
binned on ~ 30-minute intervals. We analysed a "de- trended" data 
set derived through a smoothing procedure. At first, we isolate all 
transits from the light curve. Then the moving average with a time- 



step of 0.5 days provides the mean level of the flux. Next, we con- 
struct an interpolated, reference light curve with the cubic spline on 
these nodes. Finally, we divide the raw flux, with all transits data, 
by its values of the reference, mean level flux curve. 

The de-trended data available in the MAST database exhibit 
a growth of the flux shortly before and after a particular transit. In 
some parts of the available light-curves, spanning approximately 
300 days, the measurements appear in the raw form. We did not 
use these data, aiming to analyze a possibly uniform set of obser- 
vations. 



2.1 Modeling the stellar flux 

A common model of photometric observations of a star transited 
by planetary companions consists of two major parts. The first part 
concerns the flux deficit due to small, dark objects passing in front 
of the star. At first, the average orbital periods are determined. Then 
transit depths and duration times are parametrized on the basis of 
phase-folded light-curves. Single mid-transit times are also deter- 
mined. At the second level, we can estimate the planetary masses 
and orbital elements fitting a model of motion of mutually interact- 
ing planets. 

We focus on the first level of the photometric analysis. To 
compute the flux deficit, we use the quadratic limb darkening model 
(Mandel & Agol 2002), recaUing that the Kepler- 11 light-curves 
are relatively noisy and sampled with a low frequency, 

A/(r) = 1 - yi (1 - cose) - 72 (1 - cos 9)^ , (1) 

where r is the normalized radial coordinate w.r.t. the centre of the 
stellar disk, is the angle between the direction to the observer and 
the normal to the stellar surface. The two limb-darkening coeffi- 
cients Yi and j2 must be positive and Yi + 72 < 1 (see a study of the 
limb darkening coefficients for a few target stars of the Kepler mis- 
sion, Howarth 2011). For small ratio p = Rp/Rs of planet radius 
Rp to the stellar radius Rs, Mandel & Agol (2002) found an ana- 
lytic approximation of the flux deficit, AF = AF{z; p,yi ,72) which 
depends on the normalized distance z between the centers of stellar 
and planetary disks, projected onto the sky plane (see Eq. 8 in the 
cited paper), as well as on p and yi 2- 

If more than one planet transits the star at the same time, the 
total flux deficit can be computed as the sum of the deficits caused 
by particular planets. Obviously, 71,72 are the same for all plan- 
ets, while p and z are different for each object. If transiting plan- 
ets are small, we can use a simple model of independent transits 
rather than more general treatment (e.g.. Pal 2011). Because we 
model the photometric measurements directly, by reconstructing 
the whole light-curve, we are not restricted to single transits and 
mid- transit times. Also multiple transits can be covered. In light of 
relatively narrow observational window, multiple transits are very 
helpful to constrain orbital elements of the transit model. 

Figure 1 displays a few selected fragments of the data set 
marked with red dots and error bars which are over-plotted on the 
synthetic curve best-fitting the data (blue curve). The fitting pro- 
cedure will be described in more detail in Sect. 2.3. The last panel 
shows transits of three planets (b, d and e). 

2.2 The model of orbital motion 

The orbital motion of multiple planetary system is described in 
terms of the full A/^-body problem in the Poincare reference frame 
(e.g., MorbideUi 2002). In this frame, the Cartesian coordinates 
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Figure 1. Sample synthetic light curves over-plotted on photometric measurements of Kepler-1 1. Red dots are for the observational data, blue solid curve is 
for synthetic light-curve derived in this work. The reference epoch is JD 2,455,500. 



of the planets are astrocentric, while their velocities are barycen- 
tric. The equations of motion are integrated with the second order 
symplectic integrator SABA2 (Wisdom & Holman 1991; Laskar 
& Robutel 2001). It provides 2-3 times better CPU performance 
than other algorithms, which we tested (like the Bulirsh-Stoer- 
Gragg scheme, BGS) constrained with the same time-step accu- 
racy. To speed-up the computations even more, we did not integrate 
the system at all measurements moments. This would force ^30- 
minute step- size of the integrator. Instead, we fixed this step- size to 
^ 1/20 of the innermost period of planet b, i.e., A^ ^ 0.5 day. 
Furthermore, the flux function F(t) is computed only close to 
the mid-transits. Ingress and egress times of particular events are 
tabulated. When a transit takes place, the coordinates of particu- 
lar planet at time t required to evaluate the flux deficit are deter- 
mined through the polynomial interpolation on five nodes around 
t. Through a comparison with the direct, full-accuracy integrations 
with the BGS algorithm, we found that the selected time time- step 
and the number of interpolation nodes provide a sufficient preci- 
sion and acceptable CPU overhead. We examined this method by 
changing the number of nodes in the polynomial interpolation, as 
well as the time step-size. The flux level, interpolated on five nodes 
and with A ? 0.5 day, differs from its exact value by less than 
10-9. 



2.3 Optimization algorithm and error estimation 

We searched for the best-fit model of the transits by a common 
minimization of the Xv function. This function is defined as fol- 
lows: 



1 



A^obs-A^n-1 



^obs 1 r 



-ntj) 



(2) 



where A/obs is the number of observations, A/p is the number of free 
parameters, v = A/obs — A/p — 1 is the number of the degrees of free- 
dom, Gj is the error of the j-th observation Fy, and F{tj) is a model 
function evaluated at time tj. This form of the Xv~function is cor- 
rect if the uncertainties are uncorrelated (see, e.g., Baluev 2009). 
To verify whether the available photometric data fulfill this assump- 
tion, one has to use a more general statistical model incorporat- 
ing the red-noise effect. However, under particular settings of our 
A/^-body photometric model, this would require an enormous CPU 
overhead. Hence, we use equation 2 as a reasonable first order ap- 
proximation. 



The best-fit parameters of the transits model are searched 
through a two-step optimization method. In the first step, we ap- 
ply a robust and well tested quasi-global Genetic Algorithm (GA, 
see, Charbonneau 1995; Deb 2004)^ which makes it possible to 
find promising solutions. A local, fast gradient method (here, the 
Levenberg-Marquardt algorithm) is then used to refine the solu- 
tions found in the GA step. Such an approach is called the hybrid 
optimization (see Gozdziewski et al. 2008, and references therein). 
Let us note that the parameter space is huge as it has dimension 
of 50. Some of these parameters can be determined very well, like 
the orbital periods of the transiting planets. Unfortunately, due to 
the relatively short observational time window, many parameters 
which are critical for the stability (relative inclinations, masses, 
nodal lines) cannot be well constrained. It makes the fitting pro- 
cess a challenging problem. 

The parameter errors are estimated through the bootstrap al- 
gorithm (see, e.g.. Press et al. 1992). The bootstrap is CPU- 
demanding, but it is straightforward method to estimate standard 
errors in high-dimensional problems and for large number of data. 
The light-curves which we analyzed have ^22,000 points. The 
bootstrap algorithm requires to find the best-fit solutions to a large 
number of synthetic sets derived through random sampling with re- 
placement from the original measurements. To obtain reliable error 
estimates of the best-fit parameters, one needs at least ~ 10^-10^ 
synthetic solutions. When such a large set of the best-fit models 
is gathered, we constructed normalized histograms for each free 
parameter. These histograms reflect the parameter distribution in 
response to the errors of the measurements, and may be smoothly 
approximated by an asymmetric Gaussian function. This makes it 
possible to determine the standard uncertainties. To perform the 
bootstrap procedure, at first one needs to find reliable best-fit pa- 
rameters for the nominal data set. This step was done through an 
intensive quasi-global search with the help of the hybrid algorithm. 
The bootstrap computations are CPU-time consuming and were 
performed on the reef CPU-cluster of the Poznan Supercomput- 
ing Centre. 



^ We use publicly available implementation by Kalyanmoy Deb, see http : 
/ /www . iitk . ac . in/kangal/pub . htm 
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2.4 Numerical setup of the dynamical analysis 

In spite of small eccentricities and apparently co-planar orbits, the 
Kepler- 11 systems is orbitally very active. It appears as dynam- 
ically packed planetary system (the definition is given in Barnes 
et al. 2008), with only narrow stable zones in the phase space. For 
this reason, we used the best-fit model solutions gathered in the 
bootstrap search as the input data to extensive dynamical study of 
this system. As we will discuss later, a study of the stability is a 
challenging problem. Due to relatively short observational window 
(~ 500 days), weak transits having depths comparable with the 
measurements errors and a small number of data points covering 
particular transits (typically 10 — 15), the derived initial conditions 
may be shifted away from the real configurations. 

To investigate the dynamics of the Kepler- 11 system in a 
global manner, we applied an approach in our previous papers 
which is well established in the literature. It relies on reconstruct- 
ing the structure of the phase space with the fast indicator MEGNO 
(Cincotta & Simo 2000; Cincotta et al. 2003). This dynamical char- 
acteristic makes it possible to distinguish between regular (stable) 
and irregular (chaotic, unstable) trajectories in the phase space by 
computing relatively short numerical orbits. Having representative 
solutions selected in the bootstrap statistics, we study their neigh- 
borhood on the dynamical maps. Constructing a dynamical map 
relies on two model parameters, e.g., the semi-major axes of a pair 
of planets. The selected parameters are varied in the given range at 
a discrete grid. The remaining components of the initial parameter 
vector are fixed at their nominal values. If it is necessary, they are 
altered to preserve the observational constraints. Then we calculate 
MEGNO at each point of the grid. Dynamical maps are informa- 
tive and become a standard numerical tool helpful to understand 
the global dynamics of multiple systems. 

To compute the MEGNO indicator, we must solve the varia- 
tional equations to the equations of motion of the planetary A/^-body 
problem. The Kepler- 1 1 system architecture with low-eccentric or- 
bit and small masses is an ideal target for an efficient symplectic al- 
gorithm described in (Gozdziewski 2003; Gozdziewski et al. 2008). 
The general-purpose integrators, like the Runge-Kutta or Bulirsh- 
Stoer-Gragg schemes are not efficient nor accurate enough in this 
case. These methods introduce a systematic drift of the energy and 
other integrals. To avoid such errors, and to solve the variational 
equations, we apply the tangent map introduced by Mikkola & In- 
nanen (1999). As the very basic step, it requires to differentiate 
the "drift" and "kick" maps of the standard leap-frog algorithm. 
The variations may be then propagated within the same symplectic 
scheme, as the equations of motion. Having the variational vector 
8 computed at discrete times, we find temporal y and mean Y of 
the MEGNO at the j-ih integrator step j = 1,2,..., (Cincotta et al. 
2003; Gozdziewski et al. 2008): 

y^j^ ^ (7-l)F(;-l)+);(7) ^ 

yU) = ^ja-i) + 2in(^g^) 

with initial conditions y{0) = 0, F(0) = 0, 6 = |6|. The MEGNO 
maps tend asymptotically to 

YU)=ahj + b, 

where a = 0, Z? ^ 2 for quasi-periodic orbits, a = b = for stable, 
periodic orbit, and a= {l/2)c,b = Ofov chaotic orbit with the max- 
imal Lyapunov exponent o. The tangent MEGNO map is linear, 
hence the variational vector can be normalized, if its value grows 



too large for chaotic orbits. In practice, we stop the integration if 
the MEGNO indicator reaches a given limit (usually, 7 = 5). 

The symplectic maps were propagated with the 4-th order 
SABA4 scheme in (Laskar & Robutel 2001). A choice of the 
fixed step-size must be carefully controlled. We did this, checking 
whether the relative energy error is "flat" across the dynamical map 
(Gozdziewski et al. 2008) and sufficiently small. Indeed, the step- 
size ~ 0.5 day preserved this error at a level of 10~^^ over the total 
integration times up to T ~ 40,0000 yr (~ 100,000 periods of the 
outermost planet). This time scale is long enough to detect the most 
significant 2-body and 3-body MMRs though even such integration 
period may be insufficient to detect all "dangerous" unstable reso- 
nances. Weakly chaotic motions due to multi-body MMRs still may 
lead to catastrophic events after much longer time (Gozdziewski 
et al. 2008). 

The dynamical maps in this paper have typical resolution up 
to512x512 pixels. This requires an enormous CPU-time. It is ba- 
sically not possible to perform such intensive computations on a 
single workstation. Therefore, we used our new Message Passing 
Interface (MPI) based environment Mechanic (Slonina et al. 2012) 
to perform the computations in a reasonable time in CPU-clusters^. 
They were performed on the reef cluster at the Poznan Supercom- 
puting Centre. A Mechanic run of a typical dynamical map occu- 
pied up to 1200 CPU cores for ~ 16 hours. 

2.5 Free parameters of the transit model 

The free parameters of the transit model are the stellar radius 
Rq, the limb darkening coefficients 71,72; the mass m/, radius 
Ri and orbital elements of each planet in the system, where 
/ =b,c,d,e,f,g. Planetary orbits are described through the Poincare 
geometric, osculating elements at the epoch of the first observation 
JD 2455964.51128: a tuple (a/,^/,//,^2/,C0/, f^) is for the semi- 
major axis, eccentricity, inclination to the plane of the sky, the 
longitude of ascending nodes, the argument of pericenter, and the 
mean anomaly, respectively. The orbital node of the first planet, 
i^b = 0° due to invariance of the model with respect to a rotation 
of the whole system. The inclinations are obviously close to 90°. A 
deviation from 90° is irrelevant only for single-planet systems. 

In a multi-planet system, some orbits may be inclined to the 
sky plane by angle 7^ 90°, which implies different relative incli- 
nations between orbits of particular planets, even for the same 
longitudes of nodes. Due to the invariance of transits with re- 
spect to the direction of the total angular momentum of the sys- 
tem, a combination of (// ^ 90°, ^2/) means the same geometry as 
([180° — //] ^ 90°, —Q^i). Thus, when necessary, for a given planetp 
we can fix the range of /p ^ 90° and are corrected for re- 

maining companions, in accord with the invariance relation. 

Orbital elements (a/, e/, co/, fM/) are not fully suitable for tran- 
siting systems with small relative incHnations and small eccentric- 
ities. To avoid singularities and weakly constrained elements, like 
CO/ when ei = (circular or weakly eccentric orbits), we use the 
Poincare modified elements (X/ = e/cos(0/, Yi = e/sino)/) instead 
of (^/,C0/). 

Similarly, the orbital period Pi is more suitable for model fit- 
ting than ai since the semi-major axis depends on the planetary 
mass Mi (a free parameter) and on the stellar mass mo, which is 
fixed to 0.95 m0, but it can be also fitted. Hence, we define Pi as 

^ Informations on this project may found at http : //git . astri . umk . pi/ 
pro jects /mechanic. 
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one of osculating elements related to ai through the Ilird Keplerian 
law. 

The mean anomaly Mi strongly depends on CO/. It determines 
the relative orbital phase (the mean longitude) but is also related to 
ei. This may be avoided by choosing the time of the first transit 7/ 
as a free parameter instead of Mu because it is one of the directly 
determined observables from the light-curves. Simple relations be- 
tween 7] and Poincare canonical elements may be derived easily. 

2.6 Direct and indirect transit model parameters 

The direct parameters of the transit model are determined from 
the basic observables: it is the mean period of transits , the 
depths, duration times of the transits, and shapes of the light-curves 
(through the limb-darkening coefficients). These data are usually 
derived from the period-phased light-curves of particular planets. 
The depths and durations of transit determine the ratio of planetary 
and stellar radii, Ri/Rq. If the stellar mass mo is fixed then Ri and 
Rq may be resolved. We can also determine // up to the angular mo- 
mentum direction invariance, and 7/. In general, the mean period of 
transit events P* is different from the osculating orbital period at 
the epoch of the first observation, P/. A shape of the event-period 
phased light-curves make it possible to fit the limb darkening co- 
efficients, Yi,Y2. 

These parameters of the transit model are independent on the 
the A/^-planet dynamics. Hence the remaining are indirect parame- 
ters. To resolve them, a dynamical model of the orbital evolution is 
required. The indirect parameters consist of planetary masses as 
well as orbital elements, e/, Q.i, co/ and Pi (instead of P*). Knowing 
Mi and Pi, we may fix the osculating semi-major axis ai at the date 
of the first (or prescribed) observation. 

We would like to note, that the above distinction for two 
types of model parameters is somehow arbitrary in our photomet- 
ric model. In our algorithm both the direct and indirect parameters 
are fitted simultaneously, unlike, for instance, the TTV algorithm, 
in which the direct parameters are fitted at the first stage, and the 
indirect parameters are fitted in the next step. 

Usually, the direct parameters can be estimated much more 
reliably than the indirect parameters. Even a potential derivation of 
the indirect parameters depend on the particular model of motion, 
i.e., kinematic — Keplerian, or dynamic — Newtonian, and on the 
used method of modelling the observations. In the Keplerian (kine- 
matic) model (see, e.g., Agol et al. 2005), mid-transits of a given 
planet are governed by geometric reflex motion of the star around 
the center of mass in a sub-system composed of the star and all 
inner planets. For instance, transits of planet d are affected by plan- 
ets b and c, but any outer planet does not affect transits of its inner 
companions. Hence, in accord with the Keplerian model, the indi- 
rect parameters (mg, ag, ^g, ^Ig, cOg) of the outermost planet g in the 
Kepler- 11 system cannot be determined at all. 

In a given pair of planets, the outer companion affects the tran- 
sits times of the inner planet only through gravitational mutual per- 
turbations which lead to changes of osculating orbital elements. 
To account for the mutual interactions, one has to apply the self- 
consistent A/^-body model of motion of the system. 

Usually, to resolve the indirect parameters from photometric 
observations, the well known TTV method is used (Agol et al. 
2005). It has two steps. At first, we determine the mean periods, the 
mid-transits, and then the (O-C) residua, i.e., differences between 
the measured and ephemeris transit times. The (O-C) variations are 
observables in the second step during which we search for masses, 
eccentricities, and arguments of pericenters of planetary compan- 



ions. The TTV method in this form has a limitation, because it does 
not make any use of transit depths nor their duration times. If the in- 
dividual inchnations of planets are different, the planets transit the 
parent star usually at different attitudes. Hence the transit depths 
as well as duration times may vary, like the (O-C) of mid- transits. 
This information can be used to better constrain the transit model. 

The mutual inclinations depend on the longitudes of ascending 
nodes in accord with 

cos A// J = cos// cos Ty — sin// sin Ty cos(^l/ — ^j). 

Because the inclinations of transiting planets, (//, Ij) must be close 
to 90° then Alij ^\Q.i- Q.j\. Within this approximation, the TTV 
method is apparently not sensitive for individual ^1/. In fact, differ- 
ent values of fli imply different mutual inclinations affecting the 
dynamics and (O-C). However, the dynamical variability of (O-C) 
due to mutual interactions is weaker than the geometric variability 
due to changes of transit depths and duration times reflecting the 
motion of the star around the mass center of the system. 

Overall, by direct modeling of the light-curves (photometric 
measurements), rather than the mid- transit times, we can resolve 
the (O-C) with an improved precision. Modeling the light-curves 
in terms of the A/^-body model is CPU-demanding, but it makes it 
possible to estimate individual longitudes of nodes and mutual in- 
clinations. In particular, as will be shown later, the direct method 
helped us to derive accurate relative inclinations between planets b 
and c, as well as between d and e ~ 2° ± 2°. 



3 RESULTS OF THE BOOTSTRAP ANALYSIS 

We performed the direct bootstrap TTV analysis of a few differ- 
ent orbital models of the Kepler-11 system. In the most general 
case (I), all parameters discussed in the previous section are the free 
parameters of the fit model. Some of them are poorly constrained 
by the observations, in particular, the eccentricity of planet g and 
particular longitudes of nodes. Therefore, we also studied less gen- 
eral models, in which some of weakly constrained parameters are 
fixed. In the second model (II), Xg = 0, Fg = 0, i.e., Cg = 0. In 
the third model (III), also ^Ig = 0°, while in the last model (IV), 
are all fixed at 0°. Because inclinations // are not 
exactly equal to 90°, also Alij ^ 0°. 

For each of these four transit models, we applied the boot- 
strap algorithm and we gathered sets of ~ 1500 solutions for each 
instance of the transit model. 

3.1 Model I: systems with eg^O 

Figure 2 shows an outcome of the bootstrap algorithm in the form 
of normalized histograms constructed for Xg,Fg and Q.g, and de- 
picted from the left to the right panel, respectively. The red solid 
curves illustrate the best fit asymmetric Gauss function to the his- 
togram bins. The formal lo errors are marked with red bars dis- 
played above the histograms. The best fit parameters correspond- 
ing to the maximum of the Gaussian distribution are written in the 
respective panels, and they may be compared with the nominal so- 
lutions given in Table 1. The uncertainties of the eccentricity and 
longitude of node of planet g are relatively large. 

Because the nominal system is dynamically unstable, we ex- 
amined the whole set of ^ 1500 bootstrap solutions by calculating 
their MEGNO indicator (Y) on the time interval of ~ 8000 yr. It 
corresponds to ~ 25,000 periods of the most distant companion. 
Such a characteristic time scale should be long enough to detect 
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Figure 2. Bootstrap histograms for Xg, Fg, ^Ig, transit model 1. See the text for more detail. 



unstable solutions due to low-order 2-body and 3-body mean mo- 
tion resonances (Gozdziewski et al. 2008, and references therein). 
Unfortunately, all initial configurations exhibit large values of (7), 
indicating that the system is strongly chaotic. The main source of 
instability are crossing orbits in the system, that lead to disruptive 
events , i.e., one or more of the planets were ejected from the sys- 
tem or collided with the parent star. None of the tested solutions 
passed the direct integration over 10 Myr. 

The parameter space of the Kepler- 11 system is ~ 50- 
dimensional, and the dimension of the phase space of the A/^-body 
model is 36-dimensional. The (Y) experiments indicate that this 
system can be locally chaotic and its phase space is filled with 
mostly unstable solutions. Then only small regions of stable MMRs 
may be present. In the light of a large dimension of the phase 
space, the gathered statistics of best-fit configurations is still very 
poor. We conclude that due to short data span of only ~ 500 days, 
and unconstrained elements of the most general model, we can- 
not find reliably stable solutions assuming the most general transit 
model 1. Unfortunately, in this high-dimensional problem an alter- 
nate GAMP algorithm that relies on the optimization with imposed 
stability constrains (Gozdziewski et al. 2008) would be CPU-time 
expensive. 



3.2 Transit model II: systems with eg = 

In the next model, we narrow the mostly unconstrained parameters 
of the transit model. We fix the eccentricity eg = 0, hence Xg and 
Fg are both equal to 0. The results of the bootstrap algorithm are il- 
lustrated in Figs. 3-9. All panels in these figures are constructed in 
the same manner as Fig. 2. We tested, whether the best-fit param- 
eters encompass at least marginally stable solutions with (7) ^ 2 
after T = 16000 yr. Figure 3 shows the normalized histograms for 
masses of particular planets expressed in the Earth masses. Besides 
formal uncertainties obtained through the bootstrap (filled red cir- 
cles), the best-fit parameters derived in (Lissauer et al. 2011) are 
plotted (blue filled circles). Clearly, these estimates coincide very 
well in both cases. There is one exception though, since the mass 
of planet g is not resolved in Lissauer et al. (201 1). The direct code 
helps to resolve also this mass. It is constrained surprisingly well, 
in spite of a narrow observational window. This result confirms our 
predictions. Because the orbital model is constrained by all mea- 
surements, not the TTVs only, the direct algorithm makes use of 
dynamical information contained in the transit depths and widths. 

For a reference, black and green asterisks in Fig. 3 mark 
masses of the Uranus and Neptune, respectively. The masses of 
planets b and f appear in a range specific for the super-Earths. They 
are significantly smaller than the masses of two most distant plan- 
ets in our Solar system but, as we will show in the next section, 
their chemical composition has likely much common with the ice 
giants in the Solar system. The next Fig. 4 shows histograms con- 



structed for planetary radii expressed in the unit of the Earth radius. 
These results confirm data in the discovery paper. Similarly to the 
previous plots, the radii of Uranus and Neptune are marked with 
asterisks. They are also labeled with Ry and Rn, respectively. The 
derived radius of planet g confirms a hypothesis that it may belong 
to the Uranus/Neptune-class. We note that most of the planets has 
radii smaller than Ru/n, and only planet e has its radius larger. His- 
tograms of the mean densities are presented in Fig. 5. The x-axis is 
for the density expressed w.r.t. the Earth density. Black and green 
asterisks mark the values characteristic for Uranus and Neptune, re- 
spectively. The mean densities of Saturn and water are also marked 
with the red and blue symbols, respectively. According to this plot, 
the less dense planet e has a density of Saturn. The most dense 
planet b may be almost as dense as the Earth. The densities of the 
other planets span a range characteristic for Saturn and Neptune, 
from ps to Pn- 

Figure 6 is for the bootstrap histograms constructed of the 
semi-major axes. These parameters are the best determined among 
all of the transit models, with uncertainties of the order of 10~^ au 
only. We do not compare these results with data in (Lissauer et al. 
2011) because they accounted for the formal error of the stellar 
mass. Note that we fixed mo = 0.95 mQ, because we found that this 
parameter is unconstrained by the photometric data. Yet it seems 
that the Xv('^o) function monotonically increases in the range of 
mo e (O.7,1.2)m0. 

The first five panels of Fig. 7 are for the eccentricities, and the 
bottom, right-hand panel is for AcOb,c = cOb — cOc • These histograms 
confirm that the eccentricities of planets b to f are small, typically 
less that 0.05, and the arguments of pericenters are not well con- 
strained. The last panel assures us that AcOb,c is determined with an 
error of only ^10°, recalling a narrow time- window of the photo- 
metric data. The best-fit parameters of model II are given in Tab. 2. 

Inclination 4 was constrained to the ^ 90° range, and due 
to the invariance rule implied by the direction of the total angular 
momentum, the remaining inclination // may be smaller and larger 
than 90°. We tested whether there is a correlation of the transit 
events with a given half-disc of the star. We found that both cases 
are equally possible. Because the orbits are inclined to the plane 
of the sky at angles close to 90°, the relative inclinations with the 
same longitudes of nodes may be ~ 2°-3°. As expected, the in- 
direct parameters 12/ are unconstrained, see Tab. 2. Therefore, the 
main contribution to the uncertainties of the relative inclinations 
comes from ambiguous estimates of Q.i rather than of //. 

Curiously, there appears a clear correlation between mutual 
inclinations in particular pairs of orbits, namely c and e, f and e, 
as well as d and f. This can be seen in normalized histograms 
constructed for the inclinations. Fig. 8. For a chosen planet, we 
transform // to ^ 90° range (in accord with the inclination in- 
variance rule), and we compute the bootstrap histogram for Ij. 
Panels of Fig. 8, from the left to the right, are for pairs (iJ) = 
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Table 1. Bootstrap results for transit model 1. Mass of the star is 0.95 hiq (fixed). The best-fitting stellar parameters of this model are Rq = 1.140^ 
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Figure 3. Bootstrap histograms for planetary masses, transit model II. 
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Figure 4. Bootstrap histograms for the planetary radii, transit model II. 
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Figure 5. Bootstrap histograms for the mean densities, transit model 11. 
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Figure 6. Bootstrap histograms for the semi-major axes, transit model 11. 
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Figure 7. Bootstrap histograms for eccentricities and AcL)b,c. transit model II. 
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Table 2. Bootstrap results for model II (with fixed eg = 0). Mass of the star is 0.95 hiq (fixed). Best fitted stellar parameters are Rq = 1.161^qq28, 
Yi = 0.32+[j;^^, 72 = 0.41+[j;J^, yi +72 = 0.73+[j;J^. Osculating Poincare elements are given at the epoch of the first observation JD 2455964.51128. 
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Figure 8. Bootstrap histograms for absolute inclinations, transit model II. 



(e,c),(e,f),(f,d). If /e ^ 90° then much more likely /c,/f ^ 90° 
than /c,/f ^ 90°. Similarly, if /f ^ 90°, then 1^ ^ 90° appears more 
likely than /f^ 90°. 



For particular pairs of planets, the relative inclinations can 
be determined surprisingly well. Figure 9 shows the bootstrap his- 
tograms Alij for such pairs which exhibit well constrained values. 
These histograms reveal that orbits of planets b and c are almost 
coplanar. Similarly, the pair of planets d and e form an almost 
coplanar sub-system. The mutual inclinations of orbits in these 
pairs are less than 5°, with most likely values of 2°-3°. The re- 
maining panels indicate that the mutual inclinations between five 
inner orbits remain within a few degrees range. Their upper limits 
are not so small as in the first two sub- systems. The outermost orbit 
of planet g may by highly inclined to the rest of the system, see 
errors of Q.n in Tab. 2. 



These results confirm a hypothesis in the discovery paper. In 
accord with this work, planetary orbits in the Kepler-11 system 
should be mutually inclined by no more than a few degrees. It flows 
from estimating a probability that for a given orientation of the or- 
bits, all six planets transit the star. This reasoning assumes that all 
inclinations are independent. However, we found that Kepler- 1 1 
system is composed of two or three sub-systems, which exhibit 
small mutual inclinations of orbits. Although a probability that 
the mutual inclinations between these sub- systems are significant 
seems a bit larger than for fully independent orbits, it still remains 
very small. We estimate that a randomly located observer can de- 
tect transits of all 6 planets with a probability as small as ~ 0.05%, 
for both models I and II. 



3.3 Models III (Cg = 0, = 0) and IV (eg = 0, = 0) 

The results for model III and model IV are given in Tabs. 3 and 
4. Most of these results are common for all transit models I to IV. 
There are some differences regarding a determination of the mass 
of planet g. In the realm of models III and IV (note that both have 
fixed eg = and Q.g = 0), only an upper limit of rrig smaller than 
20 — 30 Earth masses may be found. The low limits of rUg in model I 
are likely due to weakly constrained eg and Q.g. 

Let us recall that in the bootstrap set derived for model II, we 
found only two solutions with (Y) ^ 2 after 16000 yr. However, 
this integration time scale is too short to detect weak instability 
which actually leads to catastrophic disruption of these configura- 
tions. It was verified by the direct, long-term integrations. Hence, 
we did not detect any long-term stable configuration in the boot- 
strap set of model II. Similarly, stability tests performed for config- 
urations of model III did not reveal any stable models. As compared 
to model II, fixed Q.g = seem does not change the general view of 
the stability of the system. 

For model IV we found many stable configurations which con- 
firm stability analysis in Lissauer et al. (2011). They found some 
stable solutions assuming that the Kepler-11 system is strictly co- 
planar. We may conclude that a factor of small relative inclinations 
is more important for maintaining the long term stability than small 
eccentricities. This will be discussed in more detail further in this 
work. 

We examined a probability that a randomly located observer 
could detect transits of all planets in the system. This is basically 
unlikely for model III (~ 0.09%), while for model IV a probability 
of such an event is larger, and we estimate it ^ 3.4%. 
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Table 3. Bootstrap results for model III (eg = 0,^2g = 0). Mass of the star is 0.95 hiq (fixed). Fitted stellar parameters: Rq = 1.158 
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Figure 9. Bootstrap histograms for the mutual inclinations, transit model II. 



4 DISCUSSION ON THE PLANET INTERIORS 

Figure 10 shows bootstrap diagrams of a few selected pairs of pa- 
rameters. These results are for model IL The top row is for the semi- 
major axes and the planetary masses, the radii and mean densities, 
respectively. The red and green curves mark the data for Uranus and 
Neptune. The bottom row is for the mass-radius, mass-density and 
radius-density relations, respectively. Similarly, the red and green 
filled circles are for Uranus and Neptune. As we concluded above, 
the orbital solutions in set II are only marginally stable, however, it 
is a matter of unconstrained orbital angles. Note that a discussion 
in this Section concerns semi-major axes (known with an excellent 
precision) as well as planetary masses and radii. 

This figure reveals that the most inner four planets in the 
Kepler- 11 system exhibit a clear and curious anti/correlation of 
masses, radii and densities with the semi-major axes. Masses and 
radii increase with a/, while densities decrease. The last panel con- 
structed for (/?, p) shows a weak anti-correlation between the radii 
and densities, the smaller radius, the larger density. The determined 
masses and radii of the planets provide some insight into their 
chemical composition. We use a simple analytic relation between 
the radius and the mass of a cold body (Lynden-Bell & O'Dwyer 
2001; Lynden-Bell & Tout 2001) to estimate the characteristic den- 
sity po of planetary matter. This density can be compared with po 



calculated for a given number of nucleons per number of electrons 
of a chemical mixture forming the planet, jUe. The value of jUg is a 
simple mean over the elements in each chemical substance or com- 
ponent. For instance the H-He mixture has = 8/7 for the mass 
proportions 3 to 1, and ice or rock has jUe ~ 2. In this way, we can 
obtain some insight into likely chemical composition of the planets. 

Our results for Kepler- 11 are presented in Fig. 11. Black 
curves with grey areas are for po and its uncertainty Apo. Each 
panel is for one planet of the Kepler- 11 system. Data for planets b 
to planet g are displayed from the left to the right, respectively. 
The colored curves are for the Solar system, i.e., for Uranus (red), 
Neptune (green), Jupiter (blue), Saturn (violet) and the Earth (light 
blue). The density po was computed in a wide range of G [1,2]. 
These values are known relatively well for the Sun companions. 
Following Helled et al. (2011), for Uranus and Neptune one finds 
jUe ^ l.l (for the icy model) and jUe ~ 1.35 (for the rocky model). 
The density po in these particular case is plotted with filled cir- 
cles. Similarly, for Jupiter and Saturn, jje may be also estimated 
1.08 - 1.09 (Guillot 1999). Values of po for these particular jHe 
are marked with circles. Let us note that densities po of Jupiter and 
Saturn are almost identical. 

Lynden-Bell & Tout (2001) argue that po is the zero-pressure 
density po,p of the chemical mixture of the planets. Because their 
model has many simplifications, po is usually 2 — 5 times larger 
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Table 4. Bootstrap results for model IV {eg = 0, Q.i = 0, / = b, c, J, e, /, g). Mass of the star is 0.95 hiq (fixed). The best-fit stellar parameters: Rq^I. 140+qq24, 
Yi = 0.30+0 30, ^2 = 0.42+^4^, Yi +72 = OJlt^fy Osculating Poincare elements are given at the epoch of the first observation JD 2455964.51128. 
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Figure 10. The top row: mass, radius and mean density as a function of the semi-major axis (model II). Black and green solid curves are for Uranus and 
Neptune, respectively. The bottom row: mass-radius, mass-mean density and radius-mean density relations. 



then po,p. Keeping this in mind, the densities po calculated for 
Kepler- 11 planets can be compared with those of the Solar sys- 
tem planets. The value of po is best determined for planet e. Its 
is very close to the Jupiter/Saturn (J/S) value ~ 0.45 gcm~^. This 
suggests, that planet e is built mainly of a H/He mixture with mass 
proportions of the elements close to 3 / 1 with a portion of heavier 
elements contained in ices or rocks. This makes the planet classified 
as a smaller "cousin" of Jupiter and Saturn rather than of Neptune 
and Uranus, as suggested in (Lissauer et al. 2011). 

The density po of planet b is determined worse than for 
planet e. It is rather unlikely though that it belongs to the same 
class as planet e. Parameter po is larger than for Jupiter and Saturn, 
even taking into account a large uncertainty. It is also larger than po 
for Uranus and Neptune (U/N)-like planets but is smaller than po 
for the Earth. We can conclude that planet b is a small planet con- 
taining a large percentage of heavy elements in its interior, which 
is likely larger than in the ice giants. It is reasonable to classify this 
planet in the super-Earths family, although, it may be also a small 
Neptune-like planet. 

Planet f has the mass similar to planet b. However, its 
composition is likely between the Jupiter-Saturn and Uranus- 
Neptune classes. Planet d has likely similar composition as planet f. 
The best-fit estimate of po for planet c is very close to the 
Uranus/Neptune value. For planet g there is only the upper limit 



of Po- However, it is probably close to the Jupiter/Saturn value or, 
less likely, to the Neptune/Uranus value. 

These conclusions are reinforced after inspecting the bottom- 
left panel of Fig. 10 illustrating the mass-radius diagram for the 
Kepler- 11 system. The mass-radius relation computed for po and 
/Je of the Solar system planets are plotted with different colors. Data 
are shown for Uranus (red), Neptune (green), Jupiter (blue), Saturn 
(violet) and Earth (light blue), respectively. Representations of this 
relation are plotted in the mass-density diagram (the middle panel) 
and the radius-density diagram (the right panel). 

For small masses, the characteristic density po is very close to 
p (see Eq. 34 in Lynden-Bell & O'Dwyer 2001). This derived de- 
cay of p with the mean distance from the star suggests that the inner 
planets may contain larger amount of heavy elements than more 
distant companions. If this correlation can be confirmed, it may 
provide an observational constraint for the planet formation theory. 
Allowing for some speculations here, let us note that all Kepler- 11 
planets exhibit masses in the same range of a few Earth masses. 
Hence, they likely have formed in a similar way and physical envi- 
ronment (Rogers et al. 201 1). Small eccentricities and small relative 
inclinations suggest that the system evolved orbitally smoothly to- 
wards the current state, conserving the ordering of initial distances 
from the star. The observed relation between po and a/ may then 
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Figure 11. Characteristic density po of the chemical mixture of planetary interiors as functions of the mean number of nucleons per one electron, /jg (model II). 



indicate the chemical composition and mass density distribution in 
the primordial protoplanetary disk. 

We underline that the results in this section must be consid- 
ered as preliminary. Due to relatively narrow time-window of the 
photometric data, masses of the planets are determined with large 
uncertainties. 



5 RESULTS OF THE DYNAMICAL ANALYSIS 

The best-fit solutions gathered with the help of the bootstrap al- 
gorithm provide us primary information required to perform exten- 
sive study of the dynamical stability of the system. Because the or- 
bits of Kepler- 1 1 super-Earth planets are confined within the mean 
distance of Mercury in the Solar system, we could expect that the 
dynamics of this system are extremely complex. Indeed, prelim- 
inary integrations demonstrated that the Kepler- 11 system is dy- 
namically packed, in accord with a definition and the PPS hypoth- 
esis in (Barnes et al. 2008). In spite of apparently ordered config- 
urations with quasi-circular, almost coplanar orbits, and relatively 
small masses, no long-term stable model I solutions were found. 
Below, we try to resolve this paradox and we try to detect sources of 
this seemingly odd and strong instability. To illustrate the structure 
of the phase space close to the best-fit configurations, we choose 
a few representative solutions and we construct the MEGNO maps 
in their vicinity. 



5.1 Quasi-stable solutions in transit model II 

For model II, among ~ 1500 initial conditions, we found only 2 
configurations exhibiting MEGNO close to 2 after T = 16000 yr. 
Parameters of these solutions are listed in Tabs. 5 and 6. The first 
stable solution (refereed to as Ila from hereafter, see Tab. 5), has a 
relatively low mass of planet c ~ 1.5 Earth masses. Two innermost 
planets b and c have almost coplanar orbits. The next three planets, 
d, e and f also form a nearly coplanar sub-system (d-e-f), which 
is inclined to the first two orbits by large angle ^15°. The outer- 
most orbit is inclined even more, by ~ 50° to the inner subsystem 
of (b-c), and by ~ 30°to the triple-planet subsystem of (d-e-f). The 
second stable solution lib (see Tab. 6) has all masses close to the 
nominal best-fit values. The mutual inclinations of five inner or- 
bits are close to 0°, while the outermost orbit of planet g is highly 
inclined to the inner orbits by ^ 90°, similarly to solution Ila. Be- 
cause the relative inclinations between particular pairs of orbits are 
large in these best-fit solutions, such systems might be unlikely ob- 
served by a randomly located observer. We estimate a probability 



of such an event as ~ 0.07% and ^0.05% for solutions Ila and lib, 
respectively. 

5.2 Triple-planet resonances as the main source of instability 

Let us now study the vicinity of these particular solutions through 
the dynamical maps. For each initial condition of the discrete grid 
with 512 X 512 resolution, we compute (Y) over T = 8000 yr. Fig- 
ure 12 shows the MEGNO maps for solution Ila. Each panel is for 
a different pair of planets. The coordinate axes are rescaled mean 
motions centered at their nominal f?/ o- 

Xi = X 10 . 

The x/-axes span la uncertainties of the semi-major axes a/, in 
accord with Tab. 2. The semi-major axes are determined very pre- 
cisely, hence the la interval span a range of 10~^ to 10~^ au. The 
rescaled Xi are confined to order of 10. 

The MEGNO is color-coded in the dynamical maps. Blue re- 
gions mean regular solutions with (Y) ^2, while yellow color is 
for chaotic (unstable) initial conditions, (Y) ^5. The resolution is 
512 X 512 pixels, total integration time is T = 8000 yr per pixel, 
SABA4 integrator step-size is 0.5 days. A single computation of 
each map took ~ 16 hrs of 1200 CPU cores. The integrations of 
each pixel were performed up to the end time T, regardless a value 
of MEGNO. 

Still, although the maps cover tiny regions of the phase space, 
close to the fixed initial condition, they reveal a sophisticated struc- 
ture. Because we consider the dynamics in terms of conservative, 
close to integrable Hamiltonian system, this structure is governed 
by the resonant motions. A relatively short integration time ~ 10^- 
10^ characteristic periods, equivalent to the orbital period of the 
outermost planet makes it possible to detect unstable MMRs. They 
appear as yellow straight bands of different widths and slopes. 
Inspecting the dynamical maps, we can identify particular reso- 
nances. 

A condition for the mean motion resonance in the A/^-planet 
system may be written in the following form: 

Y,Pi-T7 = ^[f^^f^]^ Y.Pini= OIUJq], (3) 

i= 1 i= 1 

where rii is the mean motion of the i-ih planet, /co and fo^ are the 
fundamental frequencies associated with the pericenter arguments 
CO/ and the longitudes of nodes Q.i (for all orbits), and pi are rela- 
tively prime integers. The linear relations must obey the d' Alambert 
rule. 
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Figure 12. Symplectic MEGNO maps in the (x/,Xj)-plane for solution Ila (see the text for details). Color scale for (F) is [1,5] (blue means (F) ^ 2 and stable 
solutions; yellow is for (F) > 5 and unstable motions). Each panel is for different pair of planets labeled in its bottom-right corner. 
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Table 5. Orbital parameters of marginally stable configuration Ila. Mass of the star is 0.95 mQ. Osculating Poincare elements are given at the initial epoch 
JD 2455964.51128. 



parameter/planet 


b 


c 


d 


e 


f 


g 


m[m0] 


4.550 


1.542 


7.224 


15.698 


4.340 


18.530 


a [au] 


0.091097 


0.106515 


0.154241 


0.193939 


0.249532 


0.463815 


e 


0.00800 


0.00698 


0.00918 


0.00924 


0.01827 


(fixed) 


I [deg] 


89.048 


88.938 


89.023 


88.803 


89.339 


89.470 


^ [deg] 


(fixed) 


2.119 


15.216 


14.804 


20.626 


52.574 


CO [deg] 


158.534 


138.150 


61.831 


183.853 


296.192 


(fixed) 


^ [deg] 


47.28591 


128.23248 


118.94929 


12.27874 


151.62481 


336.22407 



Table 6. Orbital parameters of solution lib. Mass of the star is 0.95 mQ. Osculating Poincare elements are given at the initial epoch JD 2455964.51128. 



parameter/planet 


b 


c 


d 


e 


f 


g 


m[m0] 


6.227 


4.422 


8.467 


11.293 


6.866 


32.651 


a[m\ 


0.091102 


0.106525 


0.154254 


0.193942 


0.249501 


0.464000 


e 


0.02314 


0.01780 


0.01148 


0.00401 


0.01859 


(fixed) 


I [deg] 


88.000 


90.849 


89.296 


91.206 


90.677 


89.733 


^ [deg] 


(fixed) 


-1.748 


-5.944 


-2.902 


-2.393 


92.907 


CO [deg] 


178.847 


175.740 


35.793 


336.443 


350.572 


(fixed) 


M [deg] 


29.83704 


92.04062 


144.32568 


218.20059 


96.06246 


336.17121 



The two-planet MMR takes place when two values of pi are 
non-zero. If three coefficients in this linear relation are non-zero, 
it means that the system exhibits 3-body IVIIVIR. In the Kepler- 11 
system, the fundamental frequencies associated with co/ and Q.i are 
much smaller than rii (these are the secular frequencies), hence the 
right-hand sides of Eq. 3 are close to 0. This makes it possible to 
skip the secular terms, as the first order approximation, and to iden- 
tify the ]V[]V[Rs through approximate resonance conditions involv- 
ing the mean motions only. 

To identify the IVllVdRs in the 1V4EGNO maps, we apply a sim- 
ple method described in (Guzzo 2005). In the {nj^rik) -plane^, the 
slope 

dnj 

of a particular resonance line determines the ratio of coefficients pi, 
i.e., 

_Pi 
Pk 

If OLj^]^ = then the MMR forms a horizontal line, planet k is in- 
volved in the MMR, while planet j is not. If ajj^ = then the 
MMR forms a vertical line, planet j is involved in the resonance, 
while planet k is not. In these cases, other planets may be involved 
in this particular resonance. If is finite and non-zero then both 
considered planets are involved in the MMR. To identify this par- 
ticular resonance, one has to compute slopes corresponding to this 
resonance in all (n/,^^) -planes. It may be not possible, if the map 
ranges do not cover the resonance band. If a^-^^ is non-zero and fi- 
nite only for one pair of planets (j, k), it means that 2-planet MMR 
is present. It should be verified whether pjuj + p^rik ~ 0. Coef- 
ficients pj, Pk of the MMR condition can be computed from the 
slopes a^-^^. Similarly, the 3-body MMR takes place, if there exist 

^ Let us note that Fig. 12 shows the (x/,Xj) -planes but «j) may be easily 
computed from these data. 



Table 7. Three-planet resonances near the best-fit model Ila. See Fig. 12. 
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-3 
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9 
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-9 
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7 
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1* 





5 


-8 


-1 








2* 





6 


-14 


5 









relatively prime integers /, j ^ i and k ^ /, 7, such that a/j , and 
aj^k are all finite and non-zero. The 3-body MMR may be iden- 
tified by computing the slope coefficients in appropriate planes of 
the mean motions. An identification of 4-body and A^-body MMRs 
can be derived as well. 

Using this simple MMR identification algorithm, we found 
most significant MMRs close to solution Ila. The identified 3-body 
MMRs were labeled at the panels, and listed in Tab. 7. The mean 
motions of solution Ila permit a few low-order 2-body MMRs in 
the vicinity of this solution, e.g., 4nb — 5wc, l^b — 3we, In^ — Sn^, 
In^ — 2/if, 2/ie — 3/if. There is no 2-planet MMR in the range of nt 
implied by la uncertainty. All straight bands with finite and non- 
zero atj have at least two images in the planes constructed for other 
planets. All features seen in the dynamical maps correspond then 
to 3- and 4-body MMRs. Possibly, even more complex A/^-body 
resonances may be found. 

Labels in Fig. 12 corresponds to data in Tab. 7. Labels with 
asterisks are written in open circles in the dynamical maps and de- 
note MMRs in the neighborhood of solution Ila. Resonances la- 
beled with numbers without asterisks and written in filled circles in 
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Figure 13. Dynamical map of solution Ila but for the integration time 
40000 yr per pixel. Mean motion resonances are labeled in accord with 
Table 7. 

the dynamical map are also present in the neighbourhood of solu- 
tion lib. 

The MMRs {In^ - 1 IfZf + 2n^) labeled as "2" and (5wc - - 
l«e) labeled as "1*" are the most close to solution Ila. Solution Ila 
passed the 16000 yr MEGNO test as stable solution. However, be- 
cause it is close to two 3-body MMRs, and is found in a a dense 
web of low-order 3-body and 4-body MMRs, its long-term sta- 
bility cannot be guaranteed by this test. The integration time of 
16000 yr corresponds to ~ 50, 000 orbital periods of the outermost 
planet g. This time is usually too short to detect a chaotic nature of 
the orbit when the 3-body resonances are present. Unfortunately, 
the CPU-overhead does not permit to derive the dynamical map 
with the integration time per single initial condition which should 
be 10-10^-times longer. Indeed, a test run over T = 40,000 yr il- 
lustrated in Fig. 13 reveals that solution Ila is chaotic and unstable. 
Besides, almost the whole (xb,Xc)-plane corresponds to chaotic mo- 
tions with (Y) > 5. 

We analyzed the second solution lib in the same manner. The 
MEGNO maps computed for 8000 yr are presented in Fig. 14. Also 
these maps reveal a dense net of 3-body and 4-body resonances. 
Some of these resonances may be identified as close to solutions Ila 
as well as to lib. They are labeled with the same numbers written 
within filled circles, as in Fig. 12. We found also a few new MMRs, 
labeled within open circles and labeled by "3", "4" and "5". The 
remaining 4-body MMRs which are visible in this figure are not 
labeled. All identified MMRs are listed in Tab. 8. 

Similarly to configuration Ila, the integration over longer time 
of r = 40000 yr, reveals that solution lib is unstable. Almost the 
whole plane of the dynamical map corresponds to (7) > 5. The 
MEGNO map is similar to Fig. 13, and is not shown here. 

5.3 Dynamical maps in the (a)|, (Oj)- and (^i, -planes. 

Figure 15 illustrates the MEGNO maps in planes of the arguments 
of pericenters. The MEGNO is calculated over T = 8000 yr. The 
(c0/,C0j)-maps are constructed a bit differently than the mean mo- 
tions dynamical maps. Because we intent to analyse configurations 
coherent with the observations, co/ can be freely varied over the 
grid, if also the mean anomalies are modified to preserve the time 
of the first transit, Tf. For instance, for a circular orbit, when the 



Table 8. Three-planet resonances near solution lib 
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argument of pericenter is shifted from the nominal value by Aco/, 
the mean anomaly should be shifted by — Aco/. For eccentric orbits 
such a correction flows from the 1st Keplerian law. 

The MMRs form even more sophisticated web in the (co/, coy)- 
planes than in the mean motion planes. These dynamical maps re- 
veal that the stability depends on the initial arguments of pericen- 
ters. It is not obvious a priori, because eccentricities are very small. 
The initial phases of the system are preserved across the maps, 
and each point corresponds to the same 7]. Keeping in mind that 
the photometric data spanning only ~ 500 days imply weak con- 
straints on angles co/, we realize how is difficult to find a stable 
initial conditions in the huge, 50-dimensional parameter space of 
the system. Figure 16 illustrates MEGNO maps in the {ei,Q.i) - 
planes calculated over T = 8000 yr. Each panel is for one planet. 
An identification of particular MMRs is much more complex than 
in the mean motion planes. It would require the frequency analysis 
of these solutions. Still, regions of stable, quasi-periodic motions 
usually form only small islands in the phase space. For the four 
innermost planets, the regular motions are confined to only ~ 5° 
range of 12/ around the nominal value and also to a small range of 
eccentricities ~ 0.01. For the two outermost planets f and planet g 
the maps look like different. A zone of stable motions of planet g 
extends towards large Q.g. It implies a large mutual inclination of 
its orbit to the rest of the system. Small Q.g provoke unstable mo- 
tions. The nominal solution is found at the very edge between the 
regular and chaotic zone. 

5.4 Stable solution of model IV 

Let us recall that in transit model IV all = 0° . Hence, the mutual 
inclinations between all pairs of orbits are close to 0° but the system 
remains non-coplanar because // is still different from 90° (hence, 
the transits of all planets in this system can be detected by a ran- 
domly located observer with a significant probability of ~ 3.4%). 
In this case we found several solutions with MEGNO converging 
to 2 after T = 16000 yr. This indicates a possibility of quasi-regular 
orbits. We chose one of such solutions. Its parameters are listed in 
Tab. 9, and we compute dynamical maps in its vicinity. Figure 17 
shows the results of this experiment in the mean motions planes. 
The integration time is T = 8000 yr per pixel. The tested config- 
urations is found in a narrow region of regular motions. The most 
prominent dynamical feature visible in the maps is associated with 
stable 3 -planet MMR between planets b, c and d. We identified it 
as (7,-10,2,0,0,0) MMR. It has been also found in dynamical 
maps associated with solutions Ila and lib as the MMR labeled as 
"1". Due to altered parameters of the model, the unstable region of 
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Figure 14. Symplectic MEGNO maps in the (x/,Xj)-plane for solution lib. Blue color means (Y) ^ 2 (stable configuration), while yellow is for (Y) > 5 and 
unstable systems. Each panel is for different pair of planets labeled in the bottom-right corner of each particular panel. The resolution is 512 x 512 pixels, total 
integration time is 8000 yr per pixel, SABA4 integrator step-size is 0.5 days. 



Table 9. Orbital parameters of solution IVa. Mass of the star is 0.95 mQ. Osculating elements of Poincare are given at the epoch JD 2455964.51128. 



parameter/planet 



m[m0j 
a [au] 
e 

I [deg] 
CO [deg] 
[deg] 



2.359 
0.091113 
0.04423 
89.141 
20.651 



3.386 
0.106505 
0.01719 
91.215 
55.728 



5.630 
0.154243 
0.00633 
89.332 
140.753 



10.841 
0.193940 
0.00258 

88.837 
236.761 



7.524 
0.249511 
0.01073 

89.394 
355.845 



25.161 
0.463991 
(fixed) 

89.770 
(fixed) 



178.88174 209.60077 40.79259 318.51831 91.57569 336.26502 



this resonance is visible in these maps. The solid black lines plotted 
across panels shown in the top row of Fig. 17 have the slope equal to 
7/10, —7/2 and 5, from the left to the right, respectively. Other 3- 
body resonances visible in the dynamical maps form a dense web, 
which may be better seen in the (x^^Xq)- and (xci,Xf)-planes dis- 
played in the bottom panels. 

To examine the stability of this solution over longer time scale, 



we calculated a single dynamical map in the (xb,Xc) -plane for 
much longer integration time T = 40000 yr. It is shown in Fig. 18. 
This solution is located in a tiny region of stable motions. In this 
case, the 3-body MMRs has a protective role for the system, saving 
it from a disruption. Actually, this solution is chaotic. To demon- 
strate this, we computed the critical argument of the 3-body ]V[]V[R 
and we plot over two intervals of time, at the beginning of the 1 IVIyr 
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integration period (the left panel of Fig. 19) and at the end of this 
period (the right panel of Fig. 19). The critical arguments exhibit 
librations alternated with circulations. Such behavior of the critical 
argument indicates a crossing of the separatrix of the resonance, 
and chaotic dynamics. In such a case, the configuration may be 
geometrically stable over very long time but it may be suddenly 
disrupted due to a slow diffusion along the resonance. 



5.5 The Arnold web structure in the phase space 

The results of experiments described in the previous sections may 
be interpreted at the ground of the dynamical systems theory. The 
dynamical stability of planetary orbits in systems with strong per- 
turbations is influenced by even small errors and the resulting dif- 
fusion due to resonances overlapping. This dynamical phenomenon 
has been investigated in the Outer Solar system. Murray & Hol- 
man (2001) identified the chaos among the Jovian planets as result- 
ing from the overlap of the components of 3-body MMRs among 
Jupiter, Saturn, and Uranus. In spite of short Lyapunov time (10^ 



years), they found the dynamical lifetime of Uranus ^10^ years. 
In this way, the analytic theory of Murray & Holman (2001) ex- 
plained an apparent paradox of long-term stable system which is 
actually chaotic. 

The structure of the 3-body and 4-body MMRs in the Outer 
Solar system was further investigated numerically by Guzzo (2005, 
2006). Very recently, it was also studied by Quillen (2011) for 
strongly interacting extrasolar systems. Through the dynamical 
maps technique, Guzzo (2005) found that if the chaotic motions 
appear in a regular net, they may be practically stable over very 
long times. Such a state of chaotic system is called the Nekhoro- 
shev regime. On contrary, if the chaotic resonances do not consti- 
tute a regular web, and minority of orbits form a global chaotic 
zone, the stability of the system is influenced by strong chaotic 
diffusion. Such regime is related to the resonance overlap, and is 
called the Chirikov regime of the dynamics (Froeschle et al. 2000; 
Guzzo 2005). 

These results may be applied to interpret the dynamical maps 
of the Kepler- 1 1 system. The dense net of the multiple-body MMRs 
form the Arnold web in the neighborhood of the best fit model 
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Figure 16. Dynamical maps in the -plane for solution Ila. The MEGNO range [1,5] is colour coded: blue means stable solutions and yellow is for 

unstable configurations. Each map is constructed for a pair of planets labeled in the bottom-right corner. The longitudes of ascending nodes are expressed in 
degrees. The nominal solution is marked with an asterisk. See Tab. 5 for the orbital elements of this best-fit model. 




Figure 17. Dynamical maps of solution IVa in the mean motions planes. The colours and symbols are the same as in Fig. 12. 
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Figure 19. Temporal evolution of the critical argument of the 3-body MMR for the initial condition IVa. The left panel is for the beginning of the integration 
period, and the right panel is for the end of the integration period spanning 1 Myr. 




-4 -2 2 4 6 8 10 12 

Figure 18. Dynamical map of solution IVa in the mean motion plane of 
planet b and planet c but for longer integration time of T = 40, 000 yr. See 
Fig. 17 for an explanation of the symbols and colour coding. 



configurations. Our experiments reveal, that this system may un- 
dergo the Chirikov regime rather than long-term stable Nekhoro- 
shev regime. We have no strong proof of this phenomenon, be- 
cause it would require non-trivial and intensive computations of 
the chaotic diffusion. Conclusions regarding the real state of the 
Kepler- 1 1 system require more precise determination of the initial 
conditions than obtained in this work. 

We also note that solutions Ila, lib investigated in detail are 
found at the very border of the chaotic and regular zones. This 
can be well seen in the (c0/,C0j)- and (^/,^2/)-planes. A change of 
these angles constrained by the best-fit errors, could "move" the 
system into larger zones of stable motions. Because the MEGNO 
quantifies the stability of the system as a whole, such a change 
would imply a shift of the initial condition in all parameter maps 
(Gozdziewski & Maciejewski 2001). The altered initial conditions 
would be also more "distant" from the unstable strips of 3-body 
and 4-body MMRs in the mean motions planes. The dynamical 
maps would be more similar to those computed for solution IV re- 
vealing most extended zones of stable motions. Overall, the dy- 
namical state of the system is very fragile and depends on subtle 
changes of the initial conditions. 

Still most of the best-fit model configurations obtained with 
the bootstrap algorithm, which appear as regular over the short- 
term dynamics time-scale are self-disrupting. This behavior is 
similar to unstable evolution in 3-body MMRs observed in the 
HD 37124 system (Gozdziewski et al. 2008). In this sense the dy- 
namical state of the Kepler- 11 system is still puzzling. The results 



of long-term integrations of the obtained sets of initial conditions 
indicate that the system resides at the edge of dynamical stability. 
These conclusions might be changed, if more data are gathered and 
analyzed. 



6 CONCLUSIONS 

In this paper, we derived an improved method for the dynamical 
analysis of photometric light curves of stars with multiple transit- 
ing planets. Its main purpose is to determine planetary masses, as 
well as a number of indirect parameters affecting the dynamical sta- 
bility of the system. This algorithm improved the well known TTV 
algorithm. The crucial point of this method is to model the whole 
photometric curve directly with the help of an efficient symplectic 
N-ho&y integration. Such the direct approach make is possible to 
account for multiple transits, as well for the transit depths and their 
widths. This in turn makes it possible to impose dynamical con- 
straints on parameters which cannot be determined in terms of the 
TTV, like the longitude of nodes and mutual inclinations of orbits. 

With the help of this new method, we re- analyzed available 
photometric data for Kepler-11. The direct algorithm imposes con- 
straints on the mass of the outermost planet g and help us to de- 
termine the mutual inclinations between orbits of planets b and 
planet c as well as between the (d-e) pair with a good accuracy of 
2°. These results extend analysis performed in the discovery paper 
(Lissauer et al. 201 1). Overall, conclusions in this paper and in our 
work coincide very well, in spite of quite a different transit mod- 
els, optimisation algorithms and uncertainties estimation methods 
applied. 

Thanks to the in-depth analysis of the Kepler-11 light-curves, 
we investigated a possible chemical composition of the planets de- 
tected in this intriguing system. The most curious finding is a clear 
anti-correlation of the mean densities of the planets with their mean 
distances from the star. The inner planets exhibit larger abundance 
of heavy elements than the outer companions. Because all eccen- 
tricities as well as the mutual inclinations of stable systems remains 
small, the system unlikely suffered violent scattering processes in 
the past. It follows that the ordering of planets have been preserved 
since their formation. A dynamical relaxation should imply large 
Ci and Mij (see, e.g., Chatterjee et al. 2008; Adams & Laugh- 
lin 2003). These factors indicate that the primordial protoplanetary 
disk might have a significant gradient of metallicity and the present 
Kepler-11 system is a real fossil of this disk. If this suggestion is 
confirmed, it can constrain the planet formation theories, in partic- 
ular concerning multiple systems of super-Earth planets. 

This conclusion is reinforced by the dynamical analysis of the 
system. We found, in accord with the discovery paper, that the sys- 
tem is basically free from 2-planet MMRs. However, its global dy- 
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namics is governed by 3 -body and 4-body MMRs. Overlapping 
of these resonances near the best-fit solutions imply an extended 
zone of dynamical chaos and very unstable configurations. Partic- 
ular multi-body resonances may stabilize the system. We identified 
such a 3-body resonance. However, the observation window span- 
ning only ~ 500 days does not make it possible to pick up this 
MMR as the only possible. Besides, the MMRs form the structure 
of the Arnold web characteristic for the Chirikov regime. In this 
dynamical state, the phase-space orbits are strongly unstable due to 
overlapping of the MMRs. In such a case, the chaotic diffusion in 
the actions (semi-major axes) space is significant and easily leads 
to strong geometric changes of orbits. Actually, our numerical ex- 
periments indicate this in the Kepler- 11 system. It remains an open 
question though, how such an apparently ordered configuration of 
planets could survive the formation phase in extremely complex 
and fragile dynamical environment. 

The discovery team argue that the system is non-resonant. 
This factor would prevent a scenario of trapping the planets into 
MMRs during the inward migration at the early stages of the evo- 
lution. As we showed here, the system is in fact extremely resonant, 
but in quite a different sense. Combining this fact with small val- 
ues of the eccentricities and anti-correlation of the chemical com- 
position with the distance to the star, we may conclude quite an 
opposite: the migration/trapping scenario could be the only way of 
preserving the primary architecture in the present form. However, 
these suggestions might be verified only if more photometric data 
are gathered and are available. 

Most likely, many other Kepler-discovered multiple extraso- 
lar systems with transiting planets exhibit qualitatively similar be- 
haviours to that one we found in the Kepler-11. Hence, the ap- 
proach in this paper is general and may be applied in the studies 
of other compact systems composed of low-massive, super-Earth 
or Neptune-like planets. 
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